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We investigate the avalanche dynamics of the Bak-Tang-Wiesenfeld (BTW) sandpile model on scale-free 
(SF) networks, where threshold height of each node is distributed heterogeneously, given as its own degree. 
We find that the avalanche size distribution follows a power law with an exponent T. Applying the theory of 
multiplicative branching process, we obtain the exponent T and the dynamic exponent z as a function of the 
degree exponent y of SF networks as T = 7/(7— 1) and z = (y— l)/(y— 2) in the range 2 < 7 < 3 and the mean 
field values X = 1.5 and z = 2.0 for 7 > 3, with a logarithmic correction at 7= 3. The analytic solution supports 
our numerical simulation results. We also consider the case of uniform threshold, finding that the two exponents 
reduce to the mean field ones. 

PACS numbers: 89.70,+c, 89.75.-k, 05.10.-a 



Recently the emergence of a power-law degree distribution 
in complex networks, i.e., p c t(k) ~ with the degree expo- 
nent 7, have attracted many attentions ll|,|2j]. Such networks, 
called scale-free (SF) networks, are ubiquitous in nature. Due 
to the heterogeneity in degree, SF networks are vulnerable to 
attack on a few nodes with large degree |3J. However more 
severe catastrophe can occur, triggered by a small fraction of 
nodes but causing a cascade of failures of other nodes 0]. 
The 1996 blackout of power transportation network in Ore- 
gon and Canada is a typical example of such a cascading fail- 
ure 1 5 ] . As another example, malfunctioning router will au- 
tomatically prompt Internet protocols to bypass the missing 
router by sending packets to other routers. If the broken router 
carries a large amount of traffic, its absence will place a sig- 
nificant burden on its neighbors, which might bring the failure 
of the neighboring routers again, leading to a breakdown of 
the entire system eventually |6|. 

To understand such cascading failures on SF networks, we 
study in this Letter the Bak-Tang-Wiesenfeld (BTW) sand- 
pile model 1 7] as a prototypical theoretical model exhibiting 
avalanche behavior. The main feature of the model on the 
Euclidean space is the emergence of a power law with expo- 
nential cutoff in the avalanche size distribution, 

Pa(s) ~s~ z exp(-s/s c ), (1) 

where s is avalanche size and s c its characteristic size. While 
many studies of the BTW sandpile model and its related mod- 
els have been carried out on the Euclidean space, the study of 
them on complex networks has rarely been carried out. 

Bonabeau |H| have studied the BTW sandpile model on 
the Erdos-Renyi (ER) random networks and found that the 
avalanche size distribution follows a power law with the ex- 
ponent t ~ 1.5, consistent with the mean field solution in the 
Euclidean space |9]. Recently Lise and Paczuski H 1 Oil have 
studied the Olami-Feder-Christensen model 111 ill on regular 
ER networks, where degree of each node is uniform but con- 
nections are random. They found the exponent to be T « 1 .65. 
However, when degree of each node is not uniform, they 
found no criticality in the avalanche size distribution. Note 



that they assumed that the threshold of each node is uniform, 
whereas degree is not. While such a few studies have been 
performed on ER random networks, the study of the BTW 
sandpile model on SF networks has not been performed yet, 
even though there are several related applications as men- 
tioned above. 

We study the dynamics of the BTW sandpile model on SF 
networks both analytically and numerically. In the model, we 
first consider the case where threshold height of each node 
is assigned to be equal to its degree, so that threshold is not 
uniform but distributed following the power law of the de- 
gree distribution. An analytic solution for the avalanche size 
and duration distributions is obtained by applying the the- 
ory of multiplicative branching process developed by Otter 
in 1949 II 211 . The multiplicative branching process approach 
was used to obtain the mean-field solution for the BTW model 
in the Euclidean space ||9J], which is valid above the critical 
dimension d c = 4. In SF networks, due to the presence of 
nodes with large degree, the method would be useful. We 
check numerically the numbers of toppling events and distinct 
nodes participating to a given avalanches, finding that they are 
scaled in a similar fashion. Thus the avalanches tend to form 
tree structures with little loops, supporting the validity of the 
branching process approach. We obtained the exponent of the 
avalanche size distribution T = 7/(7— 1) and the dynamic ex- 
ponent z = (7— l)/(y— 2) in the range 2 < 7 < 3, while for 
7 > 3, they have mean-field values T = 3/2 and z = 2. At 
7 = 3, a logarithmic correction appears. We also performed 
numerical simulations, finding that the exponents obtained 
from numerical simulations behave similarly to the analytic 
solutions. Next we consider the case of uniform threshold 
height, obtaining that T = 3/2 and z — 2 for all 7 > 2 analyti- 
cally and numerically. 

Numerical simulations — We use the static model fl3ll to 
generate SF networks. We first start with N nodes, each of 
which is indexed by an integer i(i = l,...,N) and is assigned 
a weight equal to w, = i~ a . Here a is a control parameter 
in [0,1) and is related to the degree exponent via the rela- 
tion 7 = 1 + 1 /a for large N. Second, we select two differ- 
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FIG. 1 : The avalanche size distributions for the static model of 7 = °° 
(magenta □), 3.0 (blue A), 2.2 (green O), and 2.0 (red O). The data 
are fitted with a functional form of Eq. (1). For the fitted values of f, 
see Table I. Data are logarithmically binned. 

ent nodes i and j with probabilities equal to the normalized 
weights, Wi/Y,k w k ar, d W j/Hk w k> respectively, and attach an 
edge between them unless one exists already. This process is 
repeated until the mean degree of the network becomes 2m, 
where we use N = 10 6 and m = 2 in this work. 

Next, we perform the dynamics on the SF network follow- 
ing the rules: (i) At each time step, a grain is added at a ran- 
domly chosen node i, (ii) If the height at the node i reaches or 
exceeds a prescribed threshold z,, where we set z ( - = kj, the de- 
gree of the node i, then it becomes unstable and all the grains 
at the node topple to its adjacent nodes; 

hi^hi — ki, and hj = hj+l , (2) 

where j is a neighbor of the node i. During the transfer, there 
is a small fraction / = 10~ 4 of grains being lost, which plays 
the role of sinks without which the system becomes over- 
loaded in the end. (iii) If this toppling causes any of the ad- 
jacent nodes to be unstable, subsequent topplings follow on 
those nodes in parallel until there is no unstable node left, 
forming an avalanche, (iv) Repeat (i)-(iii). The avalanches 
without loss of any grains are regarded as "bulk" avalanches 
and taken into consideration hereafter. Note that each node 
has its own threshold, being equal to its degree, which is dif- 
ferent from the models on regular lattices. 

After a transient period, we measure the following quanti- 
ties at each avalanche event: (a) the avalanche area A, i.e., the 
number of distinct nodes participating in a given avalanche, 
(b) the avalanche size S, i.e., the number of toppling events 
in a given avalanche, (c) the number of toppled grains G in a 
given avalanche, and (d) the duration T of a given avalanche. 
To obtain the probability distribution of each quantity, we per- 
form the statistical average over at least 10 6 avalanches after 
reaching the steady state. 

The avalanches usually do not form a loop, as the probabil- 
ity distributions of the two quantities A and S behave in a sim- 
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TABLE I: Values of the avalanche size exponent T and the dynamic 
exponent z for the static model with mean degree 4 and of size 
N = 10 6 . The subscripts m and t mean the measured and theoretical 
values, respectively. Note that since the dynamic exponent diverges 
theoretically as 7 — » 2, numerical simulation data contain lots of fluc- 
tuations from sample to sample. *The case of 7 = 3 has logarithmic 
corrections in T f and Zt- 

ilar fashion. For example, the maximum area and size (A max , 
S max ) among avalanches are (5127, 5128), (12058, 12059) and 
(19692, 19692) for 7= 2.01, 3.0 and 00, respectively. So we 
shall not distinguish A and S but keep our attention mainly 
on the avalanche area distribution which has the most direct 
implication in connection with cascading failure phenomena 
in real-world networks. The avalanche area distribution fits 
well to Eq. (1), where s can represent either A or S. In order to 
check, we study the case of the ER graph, which actually is the 
case of a = of the static model, obtaining T = 1.52(1), con- 
sistent with the known result 1 8 ] . As 7 decreases from 00, % for 
7= 5.0 is more or less the same, but beyond 7= 3, it increases 
rather noticeably with decreasing 7 in Table 1. Those val- 
ues are compared with the ones obtained analytically below, 
showing a reasonable agreement. The discrepancy can be at- 
tributed to the finite-size effect. Also the probability of losing 
a grain (/ = 10~ 4 ) sets a characteristic size of the avalanche, 
roughly as s c ~ 1/ (2m/). 

It is worthwhile to note that the case of T > 1.5 has never 
been observed in the Euclidean space, suggesting that the dy- 
namics of the avalanche on SF networks differs from what is 
expected from the mean-field prediction. This feature have 
also been seen in other problems on SF networks such as the 
ferromagnetic ordering of the Ising model 1141 and the perco- 
lation problem 1 15]. 

We have also considered the avalanche duration distribu- 
tion. Since the duration of an avalanches does not run long 
enough due to the small-world effect, the duration distribu- 
tion is not well shaped numerically with finite size systems. 
Instead, we address this issue rather in an indirect manner. 
We measure the dynamic exponent z in the relation between 
avalanche size and duration, 

s ~ t z (3) 

for large t . Numerical values of z for different 7 are tabulated 
in Table I. 

Branching process — Since the quantities A and S scale in a 
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similar manner, it would be reasonable to view the avalanche 
dynamics on SF networks as a multiplicative branching pro- 
cess (kJ. To each avalanche, one can draw a corresponding 
tree structure: The node where the avalanche is triggered is 
the originator of the tree and the branches out of that node 
correspond to topplings to the neighbors of that node. As the 
avalanche proceeds, the tree grows. The number of branches 
of each node on the tree is not uniform but it is nothing but 
its own degree. The branching process ends when no fur- 
ther avalanche proceeds. In the tree structure, a daughter-node 
born at time t is located away from the originator by distance 
t along the shortest pathway. In branching process, it is as- 
sumed that branchings from different parent-nodes occur in- 
dependently. Then one can derive the statistics of avalanche 
size and lifetime analytically from the tree structure §0. 
Note that the size and the lifetime of a tree correspond to 
the avalanche size s and the avalanche duration t of a single 
avalanche, respectively. 

To be more specific, we introduce the probability that a 
certain node generates k branches, which is given by 
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I7 =1 j Pd (j)k C(y-i) 



for k> 1, 



(4) 



where is the Riemann zeta function. qo = 1 — YX=i = 
1 - C(y)/C(y- !)■ ^ Eq. a, the factor kp d (k)/£J =1 jp d (j) 
represents the normalized probability that the node gains a 
grain from one of its neighbors and 1 jk is the probability that 
the node has height k — 1 before toppling occurs. The factor 
1 jk comes from the assumption that there is no any typical 
height of a node in inactive state regardless of its degree k, 
and toppling can be triggered only when the height is k — 1 . 
The assumption was checked numerically and holds reason- 
ably well. Note that with q^ of Eq. (4), the criticality condition 
Lr=o kqk — 1 is automatically satisfied. 

Let 9>{y) = E~ = i and &{w) = U =0 ?*<»* be the 
generating functions of a tree-size distribution p{s) and q^, 
respectively. Then following the theory of multiplicative 
branching process (Hi , one finds that they are related as 



&>{(d) = G)£{& i {(d)). 



(5) 



The asymptotic behaviors of p(s) can be obtained from the 
singular behavior of J2(co) near CO = 1 of Eq. (|5j. 

The generating function £l((0) is written as £2{(d) = qo + 
Liy(o))/£(y— 1), where Liy(o) is the polylogarithm func- 
tion of order 7, defined as Liy(co) = (£»/r[y])/J°(exp(y) — 
(o)~ l y^~ l dy with the Gamma function T(y). The polyloga- 
rithm function has a branch cut [1 , °°) in the complex co-plane 
with a well-known expansion near CO = 1 11711 . As manifest 
in such branch-cut discontinuity, one can see that =S(o) is ex- 
panded near CO — 1 as 
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FIG. 2: Simulation results of the branching process with the power- 
law branching ratio with 7= 2.01 (O), 2.2 (+), 2.5 (x), 3.0 (□), 4.0 
(A), and 5.0 (0). 



to the leading order in (1 - co). Here A(y) = T(l — y)/£(y— 1) 
and B(y) = [£(/— 2)/£(y— 1)] - 1. From the relation be- 
tween £?(co) and £P(y) in Eq. (|5}, co = &(y) is obtained 
by inverting y = co/£2(co). The asymptotic behaviors of 
p(s) for large s can then be calculated through 2nip{s) = 
Jcdy^(y)y s l where C is a contour enclosing y = but not 
crossing the branch-cut [I, 00 ). We obtain that 

a{y) s-r/ir-i) (2<7<3), 

/»(.si - { bs-^Qns)- 1 / 2 (7=3), (7) 

C ( y ),-3/2 ( 7>3 ), 



where a{y) 
c(y) 



-A( 7 ) 1 /(i-r)/r[l/(l - 7)], b = y/nj6, and 
Thus the exponent T is determined 



Vl/(2*B(y)). 

to be T = 7/(7- 1) for 2 < 7 < 3 and T = 3/2 for 7 > 3. 
The value 3/2 is consistent with the mean-field value on the 
Euclidean space. This behavior of T is in reasonable agree- 
ment with that obtained by numerical simulations as tabu- 
lated in Table 1. To confirm the analytic solution for T, we 
perform numerical simulations following the branching pro- 
cess of Eq.0. Indeed, we obtain T « 2.0, 1.55 and 1.5 for 
7 = 2.01 , 3.0 and 5.0, respectively (Fig. 2). 

The distribution of duration, i.e., the lifetime of a tree 
growth, can be evaluated similarly Urill . Let r(t) be the 
probability that a branching process stops at or prior to time t. 
Then it is simple to know that r(t) = J2(r(t — 1)). For large 
t, r(t) comes close to 1 and its time variation r(t) — r(t — 1) 
is given by the right hand side of Eq. (6) with co replaced by 
r(t — 1) for each region of 7. Thus the lifetime distribution 
l(t) = r(t) — r{t — 1) is given as 



(2< 7 <3), 
(7=3), (8) 
(7>3). 
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[0,1] and when a new grain is added to a node chosen with 
probability proportional to the degree of that node. In the case 
of uniform threshold, on the other hand, we get the mean-field 
behaviors for all 7 > 2. 

The fact that T increases as 7 decreases implies the re- 
silience of the network under avalanche phenomena, by the 
role of the hubs that sustain large amount of grains thus play- 
ing the role of a reservoir. This is reminiscent of the extreme 
resilience of the network under random removal of vertices 
for 7 < 3 |3j, llSJ l2(Jl , While preparing this manuscript, we 
have learned of a recent work by Saichev et al. motivated by 
the study of earthquake avalanches |21]. Their results partly 
overlap with ours. 



FIG. 3: The avalanche size distribution for the static model with uni- 
form threshold. Shown are the case of 7 = °o (O), 3.0 (O), and 2.2 
(□). All estimated values of T are 1.50 ±0.05. 
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Since the distributions of p{s) and £(t) are originated from 
the same tree structures, p(s)ds ~ £(t)dt . Thus from Eqs. 
and we obtain the dynamic exponent defined via s ~ t z 
asz= (7- l)/(7-2) for 2 < 7 < 3 and z = 2 for 7 > 3. 
Following the same steps, we can obtain the exponents T and 
z for more general case where Zi = kf with < j3 < 1. We 
find t = (7+2/3 - 2)/(7+ j3 - 2) and z = (7+ P - 2)/(y- 2) 
for 2 < 7 < Yc ■ = j3 + 2, and T = 1.5 and z = 2.0 for 7 > 7 .. 

Sandpile with uniform threshold — We also consider the 
case that the threshold height of each node is uniform, while 
its degree is distributed following the power law. To realize 
this, we choose the threshold to be zi = 2 for vertices of de- 
gree larger than 1, and zi — 1 for those of degree 1. Then 
we modify the dynamic rule accordingly: Toppled grains are 
transferred to z, randomly selected nearest neighbors, which 
is similar to that of the Manna model 1 18]. In this case, we 
obtained T ~ 1.5 in all 7 > 2 (Fig. [3}. This can easily be un- 
derstood through the branching process analogy: In this case, 

we have qi = p d (l)/{k), qo = qi = Lk=2 ^f-\ = ^T 1 ' and 
qk = for k > 2. Thus J2((i>) is analytic for all co, yielding the 
usual mean-field exponents T = 3/2 and z — 2 for all 7 > 2. 

Summary — We have studied the Bak-Tang-Wiesenfeld 
sandpile model on scale-free networks. To account for the 
high heterogeneity of the system, the threshold height of each 
vertex is set to be the degree of the vertex. Numerical simu- 
lations suggest that for 2 < 7 < 3 the scaling behavior of the 
avalanche size distribution differs from the mean-field predic- 
tion. By mapping to the multiplicative branching process, we 
could obtain the asymptotic behaviors of the avalanche size 
and duration distributions analytically. They are described by 
novel exponents T and z different from the simple mean-field 
predictions. The result remains the same when threshold con- 
tains noise as Zi = rjjki with 77, being distributed uniformly in 



[1] R. Albert and A.-L. Barabasi, Rev. Mod. Phys. 74, 47 (2002). 
[2] S.N. Dorogovtsev and J.F.F. Mendes, Adv. Phys. 51, 1079 
(2002). 

[3] R. Albert, H. Jeong, and A.-L. Barabasi, Nature (London) 406, 
378 (2000). 

[4] D.J. Watts, Proc. Natl. Acad. Sci. USA 99, 5766 (2002). 
[5] M.L. Sachtjen, B.A. Carreras, and V.E. Lynch, Phys. Rev. E 61, 
4877 (2000). 

[6] A.E. Motter and Y.-C. Lai, Phys. Rev. E 66, 065102 (2002). 
[7] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 

(1987); Phys. Rev. A 38, 364 (1988). 
[8] E. Bonabeau, J. Phys. Soc. Japan 64, 327 (1995). 
[9] P. Alstr0m, Phys. Rev. A 38, 4905 (1988). 
[10] S. Lise and M. Paczuski, Phys. Rev. Lett. 88, 228301 (2002). 
[11] Z. Olami, H.J.S. Feder, and K. Christensen, Phys. Rev. Lett. 68, 
1244 (1992); K. Christensen and Z. Olami, Phys. Rev. A 46, 
1829 (1992). 
[12] R. Otter, Ann. Math. Stat. 20, 206 (1949). 
[13] K.-I. Goh, B. Kahng, and D. Kim, Phys. Rev. Lett. 87, 278701 

(2001) . 

[14] A. Aleksiejuk, J. A. Holyst, and D. Stauffer, Physica A 310, 260 

(2002) ; S.N. Dorogovtsev, A.V. Goltsev, and J.F.F. Mendes, 
Phys. Rev. E 66, 016104 (2002); M. Leone, A. Vazquez, A. 
Vespignani, and R. Zecchina, Eur. Phys. J. B 28, 191 (2002); G. 
Bianconi, Phys. Lett. A 303, 166 (2002). 

[15] R. Cohen, D. ben-Avraham, and S. Havlin, Phys. Rev. E 66, 
036113 (2002). 

[16] T.E. Harris, The Theory of Branching Processes (Springer- 

Verlag, Berlin, 1963). 
[17] J.E. Robinson, Phys. Rev. 83, 678 (1951). 
[18] S.S. Manna, J. Phys. A 24, L363 (1991). 
[19] R. Cohen, K. Erez, D. ben-Avraham, and S. Havlin, Phys. Rev. 

Lett. 85, 4626 (2000). 
[20] D.S. Callaway, M.E.J. Newman, S.H. Strogatz, and D.J. Watts, 

Phys. Rev. Lett. 85, 5468 (2000). 
[21] A. Saichev, A. Helmstetter, and D. Sornette, e-print 

(cond-mat/0305007l. 



